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Abstract. Attenuated Radon projections with respect to tlie weigiit function 
Wij,{x,y) = (l — x'^ — fp)f^~^^'^ are shown to be closely related to the orthogonal 
expansion in two variables with respect to W^. This leads to an algorithm 
for reconstructing two dimensional functions (images) from attenuated Radon 
projections. Similar results are established for reconstructing functions on 
the sphere from projections described by integrals over circles on the sphere, 
and for reconstructing functions on the three-dimensional ball and cylinder 
domains. 



1. Introduction 

Computer tomography (CT) offers a non-invasive method for 2D cross-sectional 
or 3D imaging of an object. In a typical CT application, the distribution of the 
attenuation coefficient through a body from measurements of x-ray transmission 
is estimated and used to reconstruct an image of the object. The mathematical 
foundation of CT is Radon transform. Let / be a function defined on the unit disk 
of the plane. A Radon transform of / is a line integral, 

(1.1) 7^e(/; t) / /(x, y)dxdy, < f? < 2^, -1 < t < 1, 

J i(e,t) 

where l{0,t) = {{x,y) : xcosd + ys'md — t} O is a. line segment inside B^. An 
essential problem in CT is to reconstruct the function / from its Radon projections. 
An algorithm amounts to an approximation to / that uses values of Tie if ] t) from 
a finite set of parameters {9,t). 

The attenuation of an x-ray beam is dependent on the energy of each photon. A 
line integral as defined in Hl.l(l represents a monochromatic x-ray. In practice, how- 
ever, an x-ray is usually polychromatic, meaning that it consists of photons with 
different energies. This could lead to artifacts in the reconstruction; see, for exam- 
ple, ^ Chapt. 4]. A polychromatic x-ray is represented by the so-called attenu- 
ated Radon projections for which the integral is taken against exp{— ae(a;, y)}dxdy, 
where ae{x,y) is a given function, instead of dxdy. Attenuated Radon transform 
appears in, for example, emission tomography jjj. The reconstruction algorithms 
for attenuated Radon data have been derived from Novikov's inversion formula f |1U| 
and |H|). See also the recent survey in 3 in this direction. 
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In the present paper we consider the special case that exp{— ae(x, y)} is given, 
or can be approximated, by the function 

(1.2) W^{x,y)^{l~x^ ^yy^-^'\ {x,y) e 

where /i > 0; in other words, as{x,y) = — (/i — 1/2) log(l — — y^). The attenuated 
Radon transform, denote by TZq, then takes the form 

(1.3) n^g{f;t):= f f{x,y)W^ix,y)dxdy, < < 27r, -1 < i < 1. 

Ji(e,t) 

Clearly this is just a special case of the attenuated Radon transform. This case, 
however, appears to be useful in understanding the effect of monochromatic and 
polychromatic x-rays. In this regard let us mention the classical example of the 
water phantom in a skull in 01 p. 121], which demonstrated that beam hardening 
causes an elevation in CT numbers for tissues close to the skull bone. The attenu- 
ated Radon transform defined in (|1.3() models the boundary behavior of the x-rays 
differently. 

Our approach is based on orthogonal polynomial expansions on . Let V^(W^) 
denote the space of orthogonal polynomials with respect to the weight function 
on B^. It is well known that 

oo oo 
k=0 fe=l 

where proj^ / is the projection of / on V^(W^). The infinite series holds in the 
sense that the sequence of the partial sums 

n 

5^(/;x,y) :=Eproj^/(a:,y), n > 0, 

converges to / as n ^ cx3 in L^{B^ ,Wfi) norm. The partial sum Snf provides a 
natural approximation to /. It turns out that there is a remarkable connection 
between S/^f and the attenuated Radon transforms, which states that 

(1.4) S2rnf{x,y) = J2 T^lif;mAt;x,y)dt, = 

^-^ I 1 ^"^ 2m + 1 

where $^ are polynomials of two variables given by explicit formulas. This repre- 
sentation provides a simple and direct access to attenuated Radon data. For the 
ordinary Radon transforms {fi — 1/2), this was discovered recently in Ap- 
plying an appropriate quadrature formula to the integrals in the expression leads 
to an approximation to / that uses discrete attenuated Radon projections. One 
important feature of the algorithm is that polynomials up to a certain degree are 
reconstructed exactly, which guarantees that the algorithm has a fast rate of conver- 
gence. Such an algorithm can be easily implemented numerically. For the ordinary 
Radon transforms, the algorithm is named OPED (Orthogonal Polynomial Expan- 
sion on the Disk) and it has proved to be a highly effective method |17[ [TH| . 

There arc other expressions in the spirit of (|1.4|l . In order to prove them, we need 
to study orthogonal expansions in terms of orthogonal polynomials with respect 
to Wf^{x,y) on B^. The case /i — 1/2 is easier since an orthonormal basis for 
V|(V7i/2) is known to be Uk{xcos +ysm -^fj), < j < k. No such convenient 
orthonormal basis is available for fj, ^ 1/2. 
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There is another advantage for considering the attenuated Radon transform 
TZg{f;t). It is known that there is a close relation between orthogonal polyno- 
mials on the unit ball and those on the unit sphere, which allows us to establish 
analogous results on the unit sphere 5^. In particular, the case = on can 
be used to show that we can reconstruct a function / from its integral projections 

(1.5) Q/(C; t)= [ /(x)dc^(x), ^ C e 5^ -1 < ^ < 1, 

where x = (xi,a;2,a;3) and duj is the surface measure on S'^. Reconstruction from 
such spherical transforms has been studied in the literature, see [HI- 

From the disk we can also extend the results to the unit ball B^ and to 
cylinder domains in M'^, taking Radon projections on parallel disks in each case. It 
turns out, however, that there is an important difference between the ball and the 
cylinder. For the cylinder domain, all results obtained in the disk can be extended 
without problem. For the unit ball, however, we still have an analogue of (|1.4|) but 
the reconstruction algorithm may no longer work as efficient as in the cylinder case. 
The problem is that the operator produced by the algorithm no longer preserves 
polynomials. 

For the algorithm on B^, we provide a numerical example in Section 2, which 
reconstructs a 2D phantom image for three different values of /i. For the transform 
on the sphere and the 3D transforms on the ball and on the cylinder domain, we 
will content with deriving the algorithms and will not discuss convergence or the 
performance of the algorithms at this time. 

The paper is organized as follows. In the following section we consider the 
reconstruction and approximation on the unit disk B^ from attenuated Radon pro- 
jections. This section is divided into several subsections, the last one includes the 
numerical example. In Section 3 the results on B^ are transplanted to those on 
the surface 5^, while the attenuated Radon projections become weighted spheri- 
cal transforms. The analogous results are then established for the unit ball B'^ in 
Section 4 and for the cylinder domain in Section 5. 

2. Reconstruction and Approximation on the unit disk 

Let n'' denote the space of polynomials of d variables and let 11^ denote the 
subspace of polynomials of total degree n in 11'', which has dimension dimllji^ = 
("d'')- ^'^^ ^ri := n^i- In this section we mainly work with the case d = 2. 

2.1. Orthogonal polynomials on the unit disk. Let Wf^ be the weight function 
defined in (|1.2|l . Let V|(W^) denote the space of orthogonal polynomials of degree 
k on B'^ with respect to the inner product 

{P, Q)iJ. = / y)Q{x, y)W^,{x, y)dxdy, a,, = {fi + l/2)/7r, 

where is the normalization constant of Wfj_, = 1/ Jga Wfj_{x)dx. Thus, P G 
VliWfj.) if P is of degree k and (P, Q)^ = for aU Q G n|_j^. We note that elements 
in a basis for V^(W^) may not be orthogonal with respect to each other according 
to our definition. A basis for V^Wfj,) is called orthonormal if the elements in the 
basis are mutually orthogonal and {P,P)^ — 1. 

The reproducing kernel of the space V^i^n) plays an important role in our 
development. In terms of an orthonormal basis {P^ : < j < A:} of V|(W^), the 
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(2.2) P,(W^^;x,y) = ^-^^-^6^_ 



reproducing kernel satisfies 

k 

(2.1) Pfe(T^^;x,y)=^P,'=(x)P,'=(y). 

i=o 

The kernel is independent of the choice of the bases of V^(VK^). In fact, a compact 
formula for this kernel can be given in terms of the Gegenbauer polynomial |13 | , 

fc + /i + 1/ 2^ 
1/2 

C^r'^' ((x,y) + v/i-||x||Vi-||ylPi) (1 - t^y'dt 

for /i > 0, the formula also holds for = upon taking limit — )■ 0. Here and 
in the following, the Gegenbauer polynomials C^(s) are orthogonal with respect to 
(l-s2)^-i/2 on [-1,1], 

(2.3) c,_i/2 j\ Cl{s)Ct{s){l - s^f-'l^ds = j^^^S,, :^ hA,, 

where C\_i/2 r(A + l)/(v^r(A + 1/2)) is the normalization constant of the 
weight function (1 — s^)^^^/^ on [—1, 1], and (a)^ :— a{a + 1) . . . (a + fc — 1). For 
fi = 1/2, Cl^~^^^^{s) — Uk{s) is the Chebyshev polynomial of the second kind. 
For the weight fmiction Wi/2{x) = 1, it is known ^ that the set 

{Uk {x cos 6 j^k +y sin Oj^k) ■ < j < k} 

forms an orthonormal basis of V^(VFi/2). The elements of this basis are the so- 
called ridge functions. In general, given an angle (j) and a polynomial p GUk :=Il\, 
a ridge polynomial is defined by 

p(0; X, y) := p{x cos (f> + j/sin 4>), cj) G [0, 27r]. 

It is easy to see that p{(l>; x, y) is a polynomial in 11^ as well. The functions 
{C^+^''^(0j- fc; X, 2/) : < J < fc}, where e,^k = J7r/(fc + 1), form a basis for V|(W^), 
abeit not an mutually orthogonal one (see, for example, 23)- The lack of or- 
thonormal ridge basis in the case of /i ^ 1/2 makes the results for attenuated 
Radon transform more difficult, as we shall see below. 

We call a polynomial P G Ilfc of one variable symmetric with respect to the 
origin if P is even when fc is even, and P is odd when fc is odd. It is known that 
C^^^^^(t) is symmetric with respect to the origin. The ridge polynomials arising 
from such a polynomial turn out to satisfy a remarkable relation. 

Proposition 2.1. For n > and k < n, the identity 

(2-4) ——J2Uk{:^;cose,sme^Pk[^;x,y)=Pki9;x,y) 

holds for all polynomials Pk G Ilfc that are symmetric with respect to the origin. 
Proof. The proof uses the following elementary trigonometric identities 



(2.5) ^sinfc^^O and ^cosk^ 



1, if fc = mod n + 1 
otherwise 
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that hold for aU nonncgative integers k. Let us prove the case k = 21. We foUow 
the proof of Proposition 2.3 in ^^1- The polynomial Pk can be written as a linear 
combination of Uk-2j for < 2j < k. Consequently, we can write P2i{S', x, y) as 

I 

(2.6) P2i{e-x,y) ^ P2i{r cos{e - 4>)) ^Y.^j{r) cos2]{e - cj,) 

3=0 

in polar coordinates x = r cos (f> and y = r sin cj), where bj (r) is a polynomial of 
degree 2j in r. Furthermore, we know that 

I 

U2i{0; COS0, sin0) ~ U2i{cos{0 — (f)) — ^ dj cos2j{0 — cj)) 

j=o 

where df) — 1 and dj = 2 for j > 1. The identities ()2.5|l and the product formula of 
the cosine function shows that 

-^^cos2j(^?-;i^)cos2j(0-;^) = <^ i cos 2j(0 -</)), ifO<z=.?<n, 



i/=0 



1, ifz = j=0. 



Let us denote by Ik the left hand side of H2.4|l . The above trigonometric identity 
implies immediately that, for < 21 < n, 

II n 

hi = E E (^) E 2z(0 - T^fr) cos m - ^) 

1=0 j=0 i/=0 

= E^j Wcos2j(^-'/') =-P2;(^cos(0 -</>)) = ^2/(e;x,y). 

This completes the proof for the case k — 21 < n. The case fc = 2Z — 1 is similar. □ 

In (|2.4II the summation is over angles, vTr/(n+l), that are equally spaced in the 
interval [0,7r). In the case that n is even, the angles can be arranged as equally 
spaced angles in [0, 2tt] by using the fact that 

l^-'J cos 2„+i - COS 2,„+i and &in 2„+i - ^m 2„+i ■ 

The result is the following proposition proved in 16 for Pk being the Chebyshev 
polynomial of the second kind. 

Proposition 2.2. For m > and k < 2m, the identity 

2m 

(2-8) T E (fe; cos 9, sin o) Pk {^,; x, y) = P^ {9; x, y) 

holds for all polynomials Pk G life that are symmetric with respect to the origin. 

There are many orthonormal bases of V^(W^) that are known explicitly (see 0). 
One that is particularly useful for us is given in terms the polar coordinates 

X = r cos (pjy = r sin <j), 0<r<l, < <j) < 27t, 

and Jacobi polynomials [21 Prop. 2.3.1]. Let pi^^'^^t) denote the orthonormal 
Jacobi polynomials, that is, 

Ca,/3 f\t^\t)pl:'f'\t)il~tril + t)^dt = S„,,„, m,n = 0,1,2,... 
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where Ca^p is the normahzcd constant so that Ca.fi ^ + t)^dt = 1. 

I,' 



Proposition 2.3. For e = or 1, define the polynomials by 



(2.9) Ptjx,y) = hJr'^-''-^'\2r^ - iy-^'Sk-2LM, 

where 

S'fe_2;,o(0) cos(fc - 200 for < 2/ < fc, 
5fe-2i,i(0) =sin(fc-2/)</> for 0<2/<fc-l, 

and 

2 r(fc - 2/ + + 3/2) 



r(/i + 3/2)r(/c-2/ + i)' 

T/ien ^/lese polynomials form an orthonormal basis for V^(Wp). 

By the definition of the reproducing kernel H2.1|l and the formula (|2.2|l . it follows 
that the above orthonormal basis satisfies 

(2.10) E PtAx,y)PUcos<t>,sm4>)^'^C^{^;x,y), 

£=04 0<2l<k 

where \ — ^ + 1/2. This formula will play an important role below. It shows, in 
particular, that the expansion of C^^^^^(0; x, y) in terms of our orthonormal basis. 
The following lemma shows the converse. 

Lemma 2.4. LetOj^u = jn/ik + l). ThenforO <2l<kife = 0and0<2l<k-l 

ife^l, 

^ c ^fl..^r^^'+l/2//)...^ ,A _ _^ ' ^ 



-Y.S,^2iA0,.k)Ci:^'^'{0j.k;^,y) = ,^ l, Ht^,d,^,PtJx,y), 

^ j=0 K + ^ -t- 2 

where di,k = 1/2 if 21 < k and di^k — ^ if2l = k, H^^, := ^P;''^^^^''^ ^''(1) '^'^^ 

(Ai + ^)/(M + i)fe-;(fc + M + i) 



/!(fc-0!(fc-^ + M+f) 



2 

Proof. Using the identities (|2.5|) it is easy to verify that 
1 ^ 

(2.11) 7— TE'5^-2^^'(^j.fc)'^fc-2;',e(^j,fe) =c^;.fc'5M'- 

Using (|2.9|1 and the fact that Pl''^{cos9j^k,sm9i^k) = Hj^i^Sk-2i,e{9j,k), we obtain 
fc 



^ E '5'fe-2;,e(6'i,fe)C^'*'^^^(6'j,fc; a;, y) 



I 1 1 



fc + /i + 



^ ^ P^fc^(a:,y)^^P,';^^(cos0,-fc,sin0,-fe)5, 



k-2Le[0j^k) 



2 0<i<2fc i=0 



/v ^ /i -j- 2 

upon using the equation H2.11(l . Finally, the expression of [i?^^.]^ is derived from 
the well-known formula of pf'^{l) (see ^Tj) and the formula of h'^j^. □ 
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Lemma 2.5. Let 9j^k be as above. Then 



1 ^ 

— — y2 Sk-2l.e{0].k)UkiOjX, cos (t),sin(l)) = di ^kSk-2l.e {((')■ 

Proof. Using (|2.t)|) and the analog formula for U21-1, the identity is an easy conse- 
quence of H2.11|l . □ 

2.2. Attenuated Radon transforms. Let 9 be an angle measured counterclock- 
wise from the positive a;-axis. Let i denote the line perpendicular to the direction 
(cos 0, sin 0) and passes through the point {t cos 6, t sin 9). The equation of the line 
is i{9,t) — {ix,y) : xcos9 + ysin9 = t} for —1 <t< 1. We use 

(2.12) I{9,t) ^ i{9,t)nB'^, 0<9<2tt, -!<<<!, 

to denote the line segment of i inside B^. Let W^. be the weight function defined 
in 1)1. 2|l . The attenuated Radon projection of a function /, with respect to W^, in 
the direction 9 with parameter t E [—1, 1] is defined in (|1.3|) . It can be written as 

(2.13) 7^^(/;^)== J ^ ^ f{tcose ~ ssm9,tsm9 + scos9)W^{s,t)ds, 

using the fact that the mapping (s,i) s- {x,y) defined by a; = tcos9 — ssm9 and 
y = isin^ -|- scos9 amounts to a rotation. When /i = 1/2, this is the usual Radon 
projection, which is also called an X-ray transform. The definition 1)1.3(1 or (|2.13(l 
shows that n^if; t) = K^gif; -t). 

The ridge polynomials are particularly useful for studying Radon transforms, as 
seen in the following result: 

Proposition 2.6. For f e L^{B'^) and p G 11^, 



(2.14) J^^ fix, y)p{^- X, y)W^,ix, y)dxdy = j ^ 7^^(/; t)p{t)dt. 

Proof. Since the change of variables t = x cos (j) + y sin and s = —x sin (f) + y cos ( 
amounts to a rotation, we have 

/ f{x,y)pk{(l)]x,y)W^{x,y)dxdy 

= / fitcoscf) ~ ss\a(l),ts\-a.(l> -\- scos(l))pk{t)Wij_{t,s)dtds 
1 /-Vi-t^ 

f{tcos(f> — s sin 0, t sin (/!> + scos (j))Wf^{t, s)dspk{t)dt, 

the inner integral is exactly TV^{f;t) by H2.13|l . □ 

In particular, attenuated Radon transforms of the orthogonal polynomials in 
V^iWfj) can be exphcitly computed. 

Lemma 2.7. If P E Vl{Wf,) then for each t e (-1, 1), < 6* < 27t, 

(2.15) 7^^(P;^) = 6^(1 - ty^^^^^P{cose,sm9), 



'k 

where ~ f^"^ '^p defined in (|2.3|l . 
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Proof. Changing variables in l|2.13(l shows that 
<3(t) :=(l-f')-"KS(f>;i) 



Since an odd power of y/l — t in the integrand is always attached with an odd 
power of s, which has integral zero, Q{t) is a polynomial of t of degree at most k. 
Furthermore, the integral shows that Q{1) — b^P{cos9,sm9). The equation (|2.14|l 
in Proposition 12 . 61 shows that 



for j — 0, 1, . . . , fc — 1, since P € Vk{B^). In particular, this shows that Q{t) is 
in fact orthogonal to all polynomials in Ilfc-i with respect to the weight function 
(1 — t^)^ on [—1, 1]. Since Q is of degree k, it must be an orthogonal polynomial 
of degree k with respect to this weight function. Hence, we conclude that Q{t) = 
cC>^+^'\t) for some constant c independent of t. Setting t = \ shows that c ~ 

b^p{cose,sme)/c^+^^^{i). a 

In the case of /i = 1/2, the above lemma appeared first in 

2.3. Orthogonal expansion and attenuated Radon projections. The stan- 
dard Hilbert space theory shows that any function in _L-^(VF^; B^) can be expanded 
as a Fourier orthogonal series in terms of V^(W^). More precisely, 

(2.16) L^{W,;B^) - £ VUW,) : / = £ projj: /, 

fc=i fc=i 

where proj^ / is the orthogonal projection of / from L?'{W^ \ B^) onto the subspace 
V^(W^). It is well known that proj^ / can be written as an integral operator in 
terms of the reproducing kernel Pk{W^ \ •, •) of Vk{B^) in L^{B^); that is, 

(2.17) proj^/(x)=/ Pfe(T4^^;x,y)/(y)T4^^(y)dy, 

where x = {xi,X2) and y = (yi,?/2)- 

This formula plays an essential role in studying the convergence behavior of the 
orthogonal expansions, see for example ilSn.lS) . For our purpose, we need a different 
expression for projj. /. This is the following remarkable formula that relates proj;, / 
to the attenuated Radon transforms of / directly. Let 

^, = —-, 0<,y<n. 
n + l 

Theorem 2.8. For n > and k < n, the operator proj^ / can be written as 

(2.18) prof, f{x,y) = a J nl{f;t)D^+'^'{^,,t;x,y)dt 

2ri+l „i 

(2.19) =^;r-^E«W nl{f;t)D'^+'^\^.,t;x,y)dt 
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where 

(2.20) Dt:'''\^,t-x^y) = ^\l\II^ Ct^''\m^''\^-.x.y) 

with Xfj^ = [Hj'J-^ and 

k 
1=0 

Proof. Since (7^^^^^ is symmetric with respect to the origin, using Proposition 
and Proposition 1^21 we have 

/ f{x,y)C'^^^''^{0jX,x,y)W^{x,y)dxdy 



1 " 

— ^Uk {^i. ; cos 0j^k , sin Bj^k) 



n , _ 

1^=0 

X / f(x,y)C'^^^^^{£,^;x,y)W^(x,y)dxdy 



— -^a^ / 7^e.(/;^)c^+'/'(^)c^^^7fe(e.;cos0, 



Using Lemma [2.41 and Lemma [2. 51 we conclude that 

/ f{x,y)Ple{^^yWt,{x,y)dxdy 
Jb^ 



= ;7^E«. / '^^M-;t)cr^\t)dt'^±^[H^^,r 

J-1 M+2 

1 

X dil I 1 'E Sk-2i.e{0],k)Uk{£.„;cos9j^k,siTi0.j,k) 
j=o 

Multiplying by Pj'^lx^y) and sum up, it follows from the definition of the repro- 
ducing kernel that 

wo]lf{x,y)=-^Y.^, / 7^,„(/;^)c^^/^(^)d^^^ 



.=0 --^ 
fc 

X 

i=0 



"^^i.=0 "'-1 '^'^s 

since P;'^^(cos^i,,sin^y) = Hif,Sk-2i.e{^v) and Af^, — [iJ^'j,]^^. This proves the first 
identity. 

We now prove the second equation p.l9|l . Using the fact that ^„+i,+i = + tt, 



cos(fc - 2l)(,n+v+i = (-1)'' cos(fc - 2l)(,^, sin(/c - 2l)£,n+u+i = (-1)'' sin(fc - 2/)^^, 



10 YUAN XU, OLEG TISCHENKO, AND CHRISTOPH HOESCHEN 

we conclude that D'^'^^^^ x,y) = (— l)''C^^^^'^(^„+i+,y; a;, y). Hence, using the 
fact that T^'^^+Tjif^i) = '^'^^if'^ conclude that 

Adding this equation and the first equation of H2.18f) and dividing the result by 2, 
we then have (EUH). □ 

In the case of /x = 1/2, it is easy to see that Aj^^, = l/(fc + 1), independent of 
Hence, for fi = 1/2, (|2.10() shows that 

and the formulas (|2.18|l and (|2.19|) are of particular simple form. This case was 
studied in 

The two expressions of projj, / look similar but are different in an important 
point. The first expression consists of Radon projections in equally spaced direc- 
tions along half of the the circumference of the circle, while the second expression 
uses Radon projections in equally spaced directions along the entire circumference 
of the circle. This distinction is meaningful for reconstruction algorithms for Radon 
data. 

If n is even, then we can use Proposition 12.21 instead of Proposition 12.11 in the 
proof. The result is another identity that uses Radon projections over equally 
spaced angles in [0, 27r]. Let 

2i/7r 

0< ly <2m. 



2m + 1' 

Theorem 2.9. For m > and k < 2m, the operator proj^ / can he written as 
(2.21) proj^/(a;,y) = -— —^a^ / Tl%^U\t)Dt'^"\4>.,t-x,v)dt 

This expression of projj. / is not a special case of H2.19II . even though both uses 
equally spaced angles. In fact, setting n — 2m shows that H2.19|l uses exactly twice 
as Radon projections in equally spaced directions. For ^ = 1/2 the identity (|2.21|l 
has appeared in The equation (|2.21|l can be deduced from (I2.18|l as follows: 
using the fact that Ti^+Trf{t) = TZ^{f; —t) and changing variable t i-^ — < in the 
integral whenever (f> = £,2v-i in H2.18|l . then making use of the equations in (|2.7|l 
and the fact that the Gegenbauer polynomial is symmetric. 

Let S!^f denote the rt-th partial sum of the expansion (|2.1t)l) : that is, 

n 

5^(/;x,2/) = Xproj^/(x,2/). 

fc=0 

The operator S"^ is a projection operator from L'^iyV^^B'^) onto H^^. An immediate 
consequence of Theorem 12.81 is the following corollary: 
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Corollary 2.10. For n > 0, the partial sum operator S^f can be written as 

1 " 

(2.22) S^{f-x,y)^—-Y,a, nl{.f-t)^^A^,,t;x,y)dt 



, 2n+l 1 



w/iere 

(2.23) <i>-(^, t; y) = X: """^t/f ^' 

fc=0 

Likewise, an immediate consequence of Theorem 12.91 is the fohowing corollary: 
Corollary 2.11. For m > 0, the partial sum operator St^^f can he written as 

(2.24) S!^^{f-X,y)^-——Y,aj 7^^_J/; i)<i>^,„ (</>,, t; a:, y)di. 

2.4. Discretization and reconstruction algorithm. The equation H2.22|l ex- 
presses the partial sum of the Fourier orthogonal expansion as the integrals of 
attenuated Radon projections in the equally spaced directions. In order to derive 
an algorithm that uses only values of attenuated Radon projections on a set of finite 
line segments, we approximate the integrals by a quadrature formula. If / is a poly- 
nomial then TV^if] t)/{l — t^Y is a polynomial of the same degree by Lemma [2. 71 
which shows that we should use a quadrature formula with respect to the weight 
function (1 — t^Y; that is, 

/ g{t){l~eYdt^Y.^,g{t,) 

where tj are real numbers and \j are chosen so that the quadrature produces exact 
values of the integrals for polynomials of degree at least M. Such a quadrature is 
said to be of N points and of precision M . A Gaussian quadrature of N points has 
the highest precision M — 2N — 1 among all quadrature formulas of N points. 

For our purpose we are interested in quadrature formulas of precision 2n that uses 
n + 1 points. A class of such formulas is given in the following proposition, which 
is based on the zeros of the quasi-orthogonal polynomial C^+y^(t) + aCn^^^'^(t), 
where a is a real number For certain range of a, such a polynomial has n + 1 
real distinct zeros in the interval [—1, 1]. 

Proposition 2.12. Let tj^n, < j < n, be the distinct zeros of a quasi- orthogonal 

polynomial C'^^^l^^{t) + aCn^^^^(t). Then there are positive numbers Aj_„ such that 
the quadrature 

(2.25) .9(0(1 - t'Ydt « hndihn) ^Hig) 

j=o 

has precision 2n if a =/= 0. If a = then the quadrature has precision 2n -|- 1. 

Using an appropriate quadrature on the integrals in H2.22() we obtain a recon- 
struction algorithm for the attenuated Radon data. We state such an algorithm 
only in the case of the quadrature formula in l)2.25|l . 
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Algorithm 2.13. Let fi > and n > 0. Let tj^n and Xj^n be as in (|2.25(l . For 

(a;, y) G define 

n n 

(2.26) ^^:(/;x,2;) = ^^7^^_^(/;^,■„)^j:,(a;,y), 

i/=0 j"=0 

where 

TU^.V) = ^^(1 -<,^)-''<fi;(e.,t,,«;a:,y). 

For a given /, the approximation process A^f uses attenuated Radon data 

{7^^^ (/; tj- „) ■.Q<v<n, < j < n} 

of /. The data consist of Radon projections on n + 1 equally spaced directions 
(specified by <^,y) along the circumference of a half circle and there are n + 1 parallel 
lines (specified by tj^n) in each direction. The algorithm produces a polynomial 
Ai^^f which is an approximation to /. In the case of /i = 1/2 the algorithm (|2.13|) 
appeared earlier in the connection to the orthogonal partial sums, however, was 
neither established nor used there. 

Theorem 2.14. The operator is a projection operator on 11^. Ln other words, 
A^J e Ul and Ai^lP ^PforP&nl. 

Proof. The function $^(^,y, x, ?/) is evidently an element in H^. It follows im- 
mediately that Al^if G n^. By definition, S/^ is a projection operator on 11^. The 
operator At^f is obtained from 5^/ by applying the quadrature (j2.25(l . exactly 
for polynomials in W^^, on (1 — t'^)~'^TZ^^{f;t)^!^{S^i,,t; •), which is a polynomial of 
degree 2n in t variable by Lemma 12.71 and 12.23|l whenever / S 11^ . Hence, the 
quadrature is exact. Thus, A'^.f = Si'J = f if f £ Ul. □ 

Alternatively, we can use a quadrature formula of proper order on the second 
expression of (|2.22(l to derive an algorithm that uses Radon projections on 2n + 2 
directions equally distributed along the circumference of the entire circle. Instead 
of stating such an algorithm we consider the case of n = 2m and use the expression 
(|2.24|l . This leads to an algorithm that sums over 2m + 1 angles that are equally 
spaced over [0, 27r], as we shall discuss in the following subsection. 

2.5. Reconstruction algorithm using attenuated Radon projections. For 

practical applications in CT, the discretization described in Algorithm 2.13 needs 
to be further specified or simplified. In fact, one has to take into consideration 
what scan geometry is used in practice. For example, the zeros of quasi orthogonal 
polynomials will not be coincide with the discrete measurement of the attenuated 
Radon projections in the usual scan geometry. If these points were used, then it 
would be necessary to introduce an interpolation process, which would introduce 
new errors. As an alternative, we suggest to use a different discretization, which 
amounts to use a different quadrature formula. 

For the ordinary Radon projections (/i = 1/2), Gaussian quadrature formulas 
for the weight function \/l — x-^ are used for the integrals in (|2.24l) to generate 
algorithms. For practical implementation in CT, the quadrature formula 

(2-27) - fit)^^ = ^yf (cos , 
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based on zeros of Tn+i{x) — cos(n + 1)9, x = cos 6, is used ^Jj. The reason for 
such a choice hes in the scanning geometry of the input data. It turns out that, for 
n = 2m, such a choice ahows us to adopt fan beam geometry and use it as parallel 
geometry in a straightforward way. 

It is possible to use the quadrature formula (|2.27() for attenuated Radon trans- 
forms t), especially when /i is a half integer. The resulted A2m will no longer 
be a projection operator, but it still reproduces polynomials of degree slightly less 
than n when is a half integer. 

Algorithm 2.15. Form > 0, {x,y) € B'^, 

2m 2m 

(2.28) At^^if; x,y) = J2Y. ^f' cos^,) Tj:,(x, y), 

where 

T^A^.y) = ■^|^^^^sinV^,$^„,(0.,cosV^j;a;,y), - ^'^12+2 ' 

The constant /i + 1/2 in Tj'^^ comes from the fact that = (/i + l/2)/7r. 
This algorithm provides an approximation for the reconstruction of a function 
f{x,y) from a set of attenuated Radon projections 

{7^^^(/;cos?^,), 0<i^<2m, 1 < j < 2m} . 

The set {(j>i, : < < 2m} consists of equally spaced angles along the circumference 
of the disk. For fi — 1/2 it has appeared in The advantage of this algorithm 
lies in the fact that it can be used with attenuated Radon data obtained from 
the fan beam geometry directly, see the discussion in jl7|. The operator, however, 
reproduces polynomials up to a lower degree. 

Theorem 2.16. Let ^ he a half integer, fi + 1/2 G N. Then the operator Af^m 
Alaorithm \2.15\ vreserves polynomials of degree 2m — 2fi; that is, ^2m-^ ~ ^ /"'^ 

Proof. The algorithm is obtained by using the Gaussian quadrature formula (|2.27|l 
to discretize the integrals in H2.24|l . that is. 



2m 

sm 



2m + 1 ^ 



in V'j^^,, (/; cos V'j )C^^^/^ (cos V'j )• 



fc=0 

If / G ^2m-2fi then using the fact that Ti-^if] t)/(l — t^Y is a polynomial of degree 
2m — 2/i, the assumption that /i is a half integer shows that 

7^0.(/;^)/yT^=(l-^'r-'/X„(/;^)/(l-^Y 

is a polynomial of 2/i — 1 + 2m — 2/i = 2m — 1. Since $2m('?!^' ^' ') i^ '^^ degree 2m 
and the quadrature 12.27|l is of precision 4m — 1, the discretization becomes exact 
in this case and we conclude that At^mf = f '^^ f ^ ^m~2n- ^ 

Let C{B'^) denote the space of continuous function on with the uniform norm 
II • llcc and let \\Al^\\ denote the operator norm of under the uniform norm. By 
A ^ B we mean that there are two constants ci and C2 such that CiA < B < C2A. 
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Evidently the convergence of the algorithm depends on ||.4^||. In fact, since in 
Algorithm 12. 131 preserves n„, it is easy to see that 

II /- <c/(i + II )£;„(/) 

where E„{f) := inf{||/— P|| : P E II^j} is the error of the best approximation of / by 
polynomials on B^. If / has r-th order continuous derivatives, then En{f) < Cfn~^ , 
in which c/ depends on the norm of the r-th derivatives of /. The same applies to 
™ Algorithm 12. 151 which preserves Ii-2m-2ti- Using the formula in H2.13|l . the 
proof of Proposition 5.1 of |16| gives the following formula of the norm of 
Algorithm EIEI 

Proposition 2.17. The operator norm ||.42,„|| of C{B^) to C{B^) is given by 

2m 2m 

\\AU\ = , max A„(x,2;), A„(x,2;) := ^ ^(sin0,- „f |Tj:,(x, y)|. 

As m ^ oo, the norm growth in an essentially polynomial order of m. Hence, 
the algorithm converges uniformly if / is sufficiently smooth. To estimate the exact 
order of Ai^^ is difficult. In the case of ^ = 1/2, it is carried out in ^B] and the 
order is ||-42m|| mlog(m+ 1). Based on this fact, we conjecture that the operator 
norm of Ai^^ is of the the order 

||yt^„|| m^+^/^log(TO + 1), asm->oo, 

which is only slightly worse than the norm ||5^|| ^ n^^^l'^ (JJJ- If the conjecture 
holds, then the algorithm will converges uniformly for smooth / S C^{B^) with 
r > fj, + 1/2. In most applications, however, the function or image could have 
jumps; that is, there is not even continuity. The numerical tests in the case of 
ordinary Radon data shows that the algorithm is stable and yields fairly accurate 
results even when the data is highly singular See also the example given in 

the following subsection. 

2.6. Numerical Example. For the numerical examples we use Algorithm 12.151 
for which the scan geometry is easy to implement. The data required are gj^i, :— 
TV^ (/; cosipj), where 4>i, = 2vk j (2m+l) stands for the 2m+l views equally spaced 
along the circumference of the region to be reconstructed and = (2j + l)/(4m+2) 
means that the x-rays in each view is distributed according to the zeros of the 
Chebyshev polynomial T2m+i- In this case the fan data can be resorted into parallel 
data (jUi). 

We reconstruct a simple analytical phantom defined by the function 

fl if 0.9 < r < 1 or < r < 0.1 

f(x,y) = < 

■'^ \q if0.1<r<0.9, 

where r = \/ + J/^, on the unit disk. This phantom contains strong singularity 
along the circles r = 0.9 and r = 0.1. The rotationally invariant nature of the 
function allows certain simplification of the algorithm. 

For the reconstruction, we choose three values of the parameter /i, /i = 0, 1/2, 3/2. 
The case /i = means the ordinary Radon transform. The case fi — means that 
the Radon transform is attenuated by the weight function (1 — — y'^)~^^^, which 
is infinity at the boundary of the disk. The case fi ^ 3/2 means that the Radon 
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transform is attenuated by the weight function 1 — a;^ — y^, which is zero at the 
boundary. In each case, the Radon data are computed analytically. 

For each of the three values of fi, we use Algorithm 12. 151 for the reconstruction 
with a moderate m = 100. The reconstructed image is evaluated on a 300 x 300 
grid. The result is shown in Figure 1 below. 




Figure 1. From left to right, ^ = 0, 1/2, 3/2. 

These images show that the function is reconstructed rather faithfully in each of 
the three cases, even though the function has strong singularity. The case fx— 1/2 
has been tested extensively and compared with FBP algorithm (^| The 
above is our first attempt to test the algorithm for attenuated Radon transforms. 

3. Reconstruction and Approximation on the unit sphere 

It is known that orthogonal polynomials on the unit ball and on the unit sphere 
are closely related ( f r2i ) . Since the approximation and the reconstruction in the 
previous section are based on orthogonal expansions on the unit disk, the relation 
suggests analogous results on the unit sphere = {{x,y,z) : + i/^ + = 1}, 
which we explore in this section. 

On the sphere we consider the attenuated spherical transform defined by 



where x = {xi,X2,xz) € S'^ , C ^ ^ and ^ 9^ 0, and dto is the measure on the 
subset {x € : (x, C) = t} which is the circle on the sphere. When /i = 0, 
this is the usual spherical transform (|1.5|) . see for example, [HI p. 33]. We will 
mainly work with the case that (^3 = 0. We say that a function is even in ^3 if 
/(a::i,a;2,a;3) = f{xi,X2,-X3)- 

Proposition 3.1. Let f be even in X3. If ( = (cos0, sin6', 0), then 

(3.1) Q^' f iC, t) =11^0 {F;t), Fixi,X2)^f(^xi,X2,^Jl-xl-xl^ 

Proof. Since / is even in X3 we have /(x) = F{xi,X2) for x G S^. The definition 
of shows that (x, Q = xi cos 9 ~[- X2 sin 6 — I{9, t). In terms of xi and X2, duj ~ 
dxidx2/ \Jl — x\ — x\. Thus, 

g^/(C;i)= / F(^x^,X2){\-x\-x'^^ '^'^^'^'^^ 



J x\ COS 6+X2 sin 6—t V "^1 "^2 

which is precisely 7?.g(i^; t). □ 



16 



YUAN XU, OLEG TISCHENKO, AND CHRISTOPH HOESCHEN 



Let Hf_i{x.) — l^al^. The space L?[H^]S'^) has an orthogonal decomposition 

oo 

(3.2) 7.2(^^.^2) ^^0^, 



fc=0 

where the subspaces Tij! contains homogeneous polynomials of degree k that are 
orthogonal to lower degree polynomials with respect to H^duj on . For /i = 0, 



nl is the space of ordinary spherical harmonics. Let 



be the orthogonal projection from L'^{H^;S^) onto 7i^. The space 7i{^ is closely 
related to the space V|(W^) discussed in the previous section (fl2 ). For our pur- 
pose, we only need the following relation on the orthogonal projections: if / is even 
in X3 then 

(3.3) Proj-H^ /(x) = projj;F(xi,a:2), 

where F is the function defined in H3.1(l . This relation, together with H3.1|l . allows 
us to express the projection operator on the sphere in terms of spherical transforms. 
Using these relations and Theorem 12. 81 we obtain the following result: 

Theorem 3.2. Let f be even in x^. For n > and k < n, the operator proj^M 
can be written as 



(3.4) 



proj^./(x) = ^ Vflp / Q''f{C;t)D^,^^^\^,,t;xi,X2)dt 



where ~ C = (cos^i,, sinfi/, 0), and D^^^"^ {^,t\x^y) is defined in ()2.20|l . 

Let YJ^f denote the n-th partial sum of the expansion l|3.2|l : that is, 

n 

^« (/;x) = ^proj„M f{xi,X2). 

The operator is a projection operator from L?{H^:S^) onto n„(S'^), the re- 
striction of on 5*^. An immediate consequence of Theorem 13.21 is the following: 

Corollary 3.3. Let f be even in X3. For n> 0, the partial sum operator Y^ ] can 
be written as 

1 " 

(3.5) i;r(/;x) = — Q''f{Cu;t)<^>';,{^,,t;xi,X2)dt 

2n+l 



1 ft-r i p 1 



where $^ is the function defined in (|2.23|l . 

For n — 2m we can also use Theorem 12.91 to derive an expression for i^2m(/)' 
which leads to the following corollary: 

Corollary 3.4. Let f be even in X3. For m > 0, the partial sum operator 1^2m/ 
can be written as 

(3.6) y^:M; x) - — — - Y^aJ Q''f{C.;tMrn{4>.,t; xi,x2)dt 



ATTENUATED RADON PROJECTIONS 



17 



where (pi, = 2m+i ' ^'^ ^ (cos 0^, sin(/i,y, 0), and is the function defined in (|2.23|l . 

In the case of /i = 0, the equations (|3.5|l and (|3.6|l are representations of the 
partial sums of ordinary spherical harmonic expansions, which are expressed in 

1 /2 

terms of the Legendre polynomial Pk{t) = {t). 

Just like the case of orthogonal expansions on the unit disk, we can use a quad- 
rature formula to obtain a reconstruction algorithm using spherical transforms. For 
example, using the quadrature formula with respect to (1 — t^)'' in ProDOsition l2.12l 
as in the case of Algorithm l2.13l we get the following result: 

Algorithm 3.5. Let f be even in X3. Let /i > 0. For n > 0, x e S'^ , 

n n 

(3.7) 5^(/;x) =^^Q^/(C.;t,,„)7;':,(a;i,a;2), 

where tj^i are as in the quadrature (|2.25|l and T^^ are defined in Alaorithm \2.1Si 
This algorithm reconstructs a function /(x) from a set of spherical transforms 
W/(C.;ij), = (cos C., sine,,, 0), Q<v<2m, 1 < j < 2m} , 

which consists of integrals over a number of circles on the sphere. These circles lie 
on planes that are parallel to the x^-asds. The circles intersect the circumference of 
a disk perpendicular to the xa-axis at equally spaced angles. The distance between 
parallel circles depends on the values of tj^n- In the case /i = 0, the algorithm 
provides an approximation to the function based on ordinary spherical transforms. 
The assumption that / is even in 2:3 implies that we can use the algorithm to 
reconstruct a function defined on the upper hemisphere from spherical transforms 
that are integrals over half circles parallel to 2:3 axis on the upper hemisphere. 

If /X is a half integer, we can also state an algorithm using the quadrature H2.27|l . 
as in Algorithm 12. 151 so that tj^n = cos jn / {2m + 1). However, in the most inter- 
esting case of /i = 0, we do not have such a somewhat simplified algorithm. 

4. Reconstruction and Approximation on the unit ball 

In this section we consider reconstruction of functions on a unit ball in 
based on the attenuated Radon projections. 

4.1. Radon projections and orthogonal polynomials. We will work with at- 
tenuated Radon projections that are integrals on line segments inside B'^ with 
respect to the weight function 

w^(x) = (1- ||x|p)^-i/2^ x = (.Ti,x2,x3) e ^i>Q. 

For our purpose, however, we will only consider those lines lying on the planes that 
are perpendicular to the X3 axis. Let 2:3 = w be such a plane. Its intersection with 
the unit ball B^ is a disk {-k: x\+ xl < Vl — w^, 2:3 = w}. A line on this disk is 
given by the equation 

^■. xcos9 + y sin 9 = t\/T~^uP, -1 < i < 1. 

Let 1(6, w; t) denote the intersection of £ with B^. The attenuated Radon projection 
on such a line is then defined by 

(4.1) K{f;t,w) := / /(x)W^^(x)dx. 

Jl{e,w;t) 
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The case /i = 1/2 again corresponds to the usual Radon projection. 

Lemma 4.1. For f g L^(Wf^; B^) and for a fixed w € [—1, 1], define a function 
on by 

9wix,y) ^ f (^Vl - w^x, \/l~w^y,u?j . 
The X-ray transform (|4.1|l is related to the 2D Radon transform by 
(4.2) n>^{f;t,w)^{l^w^rn>^{g^;t). 
Proof. Since I{9, w; t) can be represented by 

xi = -v/ 1 — w'^(t cos 9 — s sin 9) , 2:2 = Vl — ui^ {t sin 9 + t cos 9) , = w 



for s E [— Vl — t^, \n~^tP], which is a rotation around X3 axis on the plane defined 
hy X3 = w, we have 



TLe{f]t,w) = {I- w'^Y I ^gu,(icos6' - ssin6',tsin6l + scos6')W^(s,t)ds. 

The integral is precisely TZg{gw',t) by (|2.13|l . □ 

Let VniWfj,) denote the space of orthogonal polynomials with respect to Wfj, on 
B^, which contains polynomials of degree n that are orthogonal to polynomials of 
lower degrees with respect to the inner product 



{P,Q)^a^,3f P(x)Q(x)W^^(a:)dx, 

JB3 



r(M + 2) 



7r3/2r(M+ 1/2)' 



where a^^3 is the normalization constant of W^. We derive a basis for V„(W^), 
making use of an orthogonal basis for V^(W^). We note that the in these two 
notations are different, the first one is on B^ and the second one is on B^. We 
denote by Cj* the orthonormal Gegenbauer polynomial, which is equal to C^/^/h^ 
by (lOl . 



Proposition 4.2. Let {P^ : < j < k} be an orthonormal basis for Vf(VK, 



Then the polynomials 

(4.3) QLkA^,y,z) = hk{l - zY^'P^ (71=1 ' 7T=1f) ^'-'^"'■'("^ 

/or < j < fc < /, where h^ — {fi + 2)k/{fJ, + 3/2)fc, form an orthonormal basis for 

Proof. From Lemma 12.31 it is easy to see that P^ is a sum of even powers of 
homogeneous polynomials when k is even, and a sum of odd powers of homogeneous 
polynomials when k is odd. Thus, it follows that Qi^tj G . Using the fact that 
Pj is orthonormal, it follows from the integral relation 



(4.4) / /(x)dx= /' / 



/ ( xi^/l - xl,X2\ll - xl,X3 ] dxidx2{l - x^)dx3 

52 
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and the fact that a^,3 — where is the normahzation of Wfj, on and 

is defined in (|2.3|l . that 



M,3 / (9/.fc.j(x)(3/',fc'j'(x)T4^^(x)dx 

JS3 



'-1 

,2 '^M + 1/2 r r jr 

— /Ife OlJ'Ok^k'OjJ' ■ 

Cfc+p+1/2 

It follows from the definition of that c^+i/2/cfe+p+i/2 = (/i + 3/2)k/ip + 2)fe, 
which completes the proof. □ 

The attenuated Radon transforms of this basis can be computed explicitly. 
Proposition 4.3. Let n > and let Qi.k.j he defined by (|4.3|l . Then 

^ ■ ' {i~t^y{i-w^Y 

= f; , -I <3/,fcj [V^-w'^ COS0, \/l - u;2 sine/), w) . 

(1) 

Proof. By Lemma l4.1l and the definition of Qi^kj we have 

where gw{x,y) = /ifcP/(x, y)(l - ■u;2)''V2c'^fc+A'+i(u;), By Lemma it follows that 

7^^(ff„;^) = - w^f'^ct+i:+\w)n>;{p^-,t) 

= ^(1 - i'T ^+1/2 Q^fcJ (Vl-w2 COS0, Vl^^sin0, w;j 

by the definition of Qi.k.j ■ Putting these equations together completes the proof. □ 

Let projfg denote the projection operator from L^{W^\B^) onto the space 
Vf{W^). Again we have the decomposition 

(4.6) L\W,-B^) = £ VliW,) : / = wofk,, /■ 

fc=0 fc=0 

Proposition 4.4. For n > and < I < n, 



(4.7) proji;3/(x) = -^^ J J ^TZl{f;t,w)Gi{^,,t,w;^)dtVT 



— w^dw 

T9 £ / / T^tif-^t'^)Gi{^-'t^^-'^)dtVT^dw 



iy=0 ' 
2ri+l „i „i 



2n + 2 



-1 J-1 
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where 

I 



Proof. The projection operator has an integral expression just as that of l|2.17|l . 
Furthermore, the kernel function P(Wp;x, y) can be written as a sum of an or- 
thonormal basis. In particular, 

I k 

Proji^a / W = X! X! fi,k,3Qi,k,] (x), 
where Qi^k,j is the orthonormal basis for Vf{W^) defined in 1)4.3(1 and 



fi,kj^af,^3 / /(y)Q(,fe,j(y)Wp(y)dy. 
Using (|4.4|l . the definition of Qk.ij, and the fact that — we have 

= c^+1/2 / / gw{u,v)Pj{u,v)W^{u,v)dudv 

where g^, is defined as in Lemma ITTI Hence, it follows from 1(2.17)1 and ((2.1(1 that 
proj,,3 /(x) = J2 hlCt^i:+\x,){l - xlf 1^0,^,1^ 

k=Q 

The identity 1(4.7(1 follows from the above equation upon using ((2.18(1 and 1(4.2(1 . □ 
Let us denote by Sl^^f the n-th partial sum of the orthogonal expansion ((4.6(1 . 

n 

5^,3/(x)=EpK3/W- 

As an immediate consequence of Proposition 14 . 41 we have 
Corollary 4.5. For n> 0, 

(4.8) 5^3/(x)-— / 7^^^(/;^,^«)ci>(^(C.,t,^«;x)d^y^^d^« 

^ 2n+l „i „i 

where 



$^i(e,i,u;;x) =^Gi(e,t,zi;;x). 



In the case of n = 2m we can use ((2.21(1 instead of 1(2.18(1 in the last step of 
the proof of Proposition 14.41 to get an expression for projj'g /. The corresponding 
expression for the partial sum is the following result: 
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Proposition 4.6. For m > 0, 

+ ^^qJ-iJ-i 

From such an expression of S'^ 3 we naturally want to derive an algorithm as in 
the 2D case. However, there is a problem when we use quadrature formula. Indeed, 
in order to obtain an algorithm, we need to discretize the integrals 

(4.9) / 1 / 1 ^^-^ ^' - w^dw 

in Sl^ 3/ by a quadrature formula. We can use, for example, the quadrature (|2.25|) 
of precision 2n, which we denote by 

/I " 
f{t){l-t'rdt^Y.^lnf{tln) 

"1 k=0 

to emphasis the dependence of tk^n and Afc.„ on the weight function. If we follow 
the 2D case, then the equation (|4.5|) indicates that we should apply the quadrature 
with respect to (1 — i^)^ in t variable, and apply the quadrature with respect to 
(1 — ti;2)'^+i/2 in w variable. The result of using these quadrature formulas gives 
the following: 

Algorithm 4.7. Let^i>0. For n > 0, x — {xi, X2,X3) e , 

n n n 
u=Oj=0 k=0 

where 

However, this is likely not an accurate algorithm. The problem is that the 
operator does not preserve polynomials of degree n. In fact, in order that 
BnP = P for P g H^, we need the discretization of the integrals l|4.9|l to be exact 
whenever / is a polynomial of degree at most n. The function 

F^(i, w) (1 - _ w^Y'^nl (/; t, u;)$^J(C., t, w; x) 

is a polynomial of degree 2n in variable t whenever / is a polynomial of degree n 
by the definition of $^ and Proposition ^31 so that the discretization in t variable 
is exact. However, the function F^{t, w) is not a polynomial in w variable. By the 
definition of Qi.k.j in the equation H4.5|l shows that F^{t,w) with / = Qi^k,j 

contains {I ~ wf/'^Ci^j^"^^ {w), which is not a polynomial in w variable if k is odd. 
The formula of $J^(Ci/, w; x) shows that it is a sum of functions, which is also not 
a polynomial. This means that the quadrature will not be exact and polynomials 
are not preserved by S^ . 

An algorithm should have high convergence order if it preserves polynomials up 
to certain degrees. The fact that B!^f does not preserves polynomials means that 
the convergence of the algorithm may not be as desirable. 
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5. Reconstruction and Approximation on the cylinder domain 

In contrast to the unit ball in M'^, the reconstruction algorithm on a cylinder 
domain works well. Let i > and let be the cylinder domain defined by 

Bl = B^x [0,L]={{x,y,z) : {x,y) G B^,0<z<L}. 

We will show that the partial sum operator of the orthogonal expansions on Wl 
admits an expression that relates to Radon data and use it to get a reconstruction 
algorithm. 

Let be defined as in (|1.2II . Let W^i.l be the weight function 

W^,.L{x,y,z) = W^,{x,y)WL{z), {x,y,z) G Bl- 

We retain the notation Til^{g\t) for the attenuated Radon projection of a function 
(7 : I— !■ R, as defined in (|1.3|l . For a fixed z in [0, L], we define 

(5.1) 7^^(/(•,•,z);^) / f{x,y,z)W^{x,y)dxdy, 

which is the attenuated Radon projection of / in a disk that is perpendicular to 
the z-axis. 

We consider the orthogonal polynomials with respect to the inner product 
(5-2) f{x,y,z)g{x,y,z)W^,L{x,y,z)dxdydz. 

Let V^(1V^_l) denote the subspace of orthogonal polynomials of degree n on Bj^ 
with respect to the inner product (|5.2|) : that is, P G Vi^(VF^,L) if {P,Q)bl — for 
all polynomial Q £ ^n^i- 

Let pk be the orthonormal polynomial of degree n with respect to Wl on [0, L] 
and let {Pj'{x, y) < j < k} denote an orthonormal basis of Vl{Wf_i). Since W^,l 
is a product on a product domain, the following proposition is obvious. 

Proposition 5.1. An orthonormal basis for Vf{Wfj_^L) is given by 

= {Plu^, : < j < fc < n} , V, z) = Pf{x, y)pn^k{z). 

In particular, the set {P; : < I < n} is an orthonormal basis for H^^ . 

For / e L^{W^,L] Bl), the Fourier coefficients of / with respect to the orthonor- 
mal system {P; : / > 0} are given by 

flM,, = «M / fi^)PtM,, W'^x. < J < fc < /. 

Let S^^ j^f denote the Fourier partial sum operator, 

71 I k 
1=0 k=0 3=0 

Just like its counterpart in two variables, this is a projection operator. The following 
is an analogue of Theorem 12 . 1 01 for the cylinder domain Bl- 

Theorem 5.2. For n> 0, 



(5.3) 5,':.^/(x) = ^Ea^ r f nl{f{;;w);t)^^^{i., 



w, t; y.)WL {w)dw dt 
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where 

(5.4) cl>jj(e, t; x) = 5] +1/2 E P^MpK^s)- 

fe=0 I 1^0 



Proof. By the definition of fj^f, j we can write 



ftkj^(^t^ / /i- fe (a; , y ) -P,- (a; , y ) H/^ (a; , y)dxdy 



where 

fL 



fi-kix,y):= / f{x,y,w)pi-k{w)WLiw)d'w, l>k>0. 
Jo 

Consequently, by the definition of proj^! in (|2.17|l , it follows that 



n I 



1=0 k=0 

We can then use the expression (|2.18() for proj^ / and the fact that 



/ n^{f{-,-,w);t)pi^k{w)WL{w)dw 
Jo 

to complete the proof. □ 

In the case of n = 2m, we can use H2.21II in place of H2.18|l in the proof. The 
result is the following proposition which has appeared in ^16j when fi = 1/2. 



Proposition 5.3. For m >0, 



(5-5) S^nuLfi^) 



, 2m „i „L 

,,_n J —1 JO 



w, t; x) Wi {w)dw dt. 



From the expression (|5.3|) or H5.5|) of 5,^^ we can apply a quadrature formula 
to get a reconstruction algorithm on Bi for the attenuated Radon data. In J2| the 
weight function Wl is chosen to be the Chebyshev weight function 

Wl{z) = -^^=, ^e[0,L], 

normalized to have integral 1 on [0, i]. The reason for this choice is that the 
Gaussian quadrature formula takes a simple form 

1 " 1 ^ 

(5.6) / g{z)WL{z)dz^—-Y,g{zi), = (l + cos 

which is of precision 2n + 1. We can apply this quadrature for the integral with 
respect to w and use the quadrature H2.25|l for the integral with respect to t in (|5.3|) 
or 1)5. 5|) . The result is the following algorithm: 

Algorithm 5.4. Letfi>0 and let jf^j^i — TZ'^ (/(•,•, z^); fj,„). Forn>0 

(5.7) s,':.l(/;x) = EEE7..,.r.^,.:(x) 

i/=0 j=0 i=0 
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where 

Like the algorithms in the previous sections, this algorithm produces a polyno- 
mial as an approximation to the function. It does preserve polynomials of lower 
degrees. 

Theorem 5.5. The operator B'^ ^ is a projection operator on 11^. In other words, 
Bnf e Ul and BnAf) = f ^ff&^l- 

Proof. Let P^ ^, j be defined as in Proposition 15.11 It follows from the definition 
in (EH that 'R-'l{Pt,k,ji-^ ■,w);t) = 'R^^{Pf\t)pi_k[w). Consequently, it follows from 
(I2.15|l that 7?.^(P(-, •, w); i)/(l — t^Y is a polynomial of degree n in both t vari- 
able and w variable whenever P G 11^. By its definition in (|5.4() . the function 
(^, w, t; x) is evidently a polynomial of degree n in both t and w variables. Hence, 
we can apply H5.6|) for w variable and apply the quadrature (|2.25|) of precision 2n 
to t variable, which are exact on (1 - t'^)-f'n'^{P{-, •, w); w, t; •). □ 

The approximation process in Algorithm 15.41 uses the attenuated Radon data 

{7^^_^(/(•,•,z,);^J■„) :0<i^<n, < j < n, < i < n} , 

which consists of Radon projections on n + 1 disks that are parallel to the z-axis. 
In other words, it consists of reconstructions of the function on n + 1 planes. 

In the case of n = 2m and /i is an half integer, we can also use the quadrature 
()2.27|l to derive a more explicit algorithm as in Algorithm 12. 151 Such an algorithm 
is given in (I6j for /z = 1/2. We shall not elaborate further. 
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